In [16]:
%pylab inline
import pydisdrometer
In [17]:
filename = '../testdata/IFloodS_APU01_2013_0502_dsd.txt' #Parsivel 05, March 2nd
In [18]:
dsd = pydisdrometer.read_parsivel_nasa_gv(filename)
So at this point we have the drop size distribution read in. NASA strips out rainfall information though and I have not yet re-added the calculation for that(Maybe today or tomorrow). Let's do some T-Matrix scattering though. This should take a little bit(up to a minute or so on my laptop depending on the file).
In [19]:
dsd.calculate_radar_parameters() #Ignore the warnings, we don't do QC first so sometimes doing scattering on zero drops causes complaints.
This assumes BC shape relationship, X band, 10C. You can pass in a new wavelength to change that. Let's plot some of the parameters, and then try to do something with the data.
In [20]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.plot(dsd.time/60.0, dsd.fields['Zh']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Reflectivity(dBZ)')
plt.xlim(5,24)
plt.subplot(2,2,2)
plt.plot(dsd.time/60.0, dsd.fields['Zdr']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Differential Reflectivity(dB)')
plt.xlim(5,24)
plt.subplot(2,2,3)
plt.plot(dsd.time/60.0, dsd.fields['Kdp']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Specific Differential Phase(deg/km)')
plt.xlim(5,24)
plt.subplot(2,2,4)
plt.pcolor(dsd.time/60.0, dsd.diameter, np.log10(dsd.Nd.T), vmin=0, vmax=6)
plt.xlabel('Time(hrs)')
plt.ylabel('Diameter')
plt.colorbar()
plt.ylim(0,5) #Zoom in some
plt.xlim(5,24)
plt.show()
Let's take what might be another common use case. Finding areas with certain size drops. In this case everything above 4mm. There are a few ways to do this but let's just create an indicator function for now.
In [21]:
plot(dsd.time/60.0 , 0<sum(dsd.Nd[:,dsd.diameter>4], axis=1))
plt.xlim(5,24)
Out[21]:
Or maybe we just want a listing of the time steps where there are drops greater than 4mm
In [22]:
time_hrs = dsd.time/60.0
print(time_hrs[0<sum(dsd.Nd[:,dsd.diameter>4], axis=1)])
In [22]: